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ON THE UNIVERSALITY OF THE KOLMOGOROV CONSTANT IN NUMERICAL 

SIMULATIONS OF TURBULENCE* 

P.K.YEUNGt AND YE ZHOU* 


Abstract. Motivated by a recent survey of experimental data [K.R. Sreenivasan, Phys. Fluids 7, 
2778 (1995)], we examine data on the Kolmogorov spectrum constant in numerical simulations of isotropic 
turbulence, using results both from previous studies and from new direct numerical simulations over a range 
of Reynolds numbers (up to 240 on the Taylor scale) at grid resolutions up to 512 3 . It is noted that 
in addition to scaling, identification of a true inertial range requires spectral isotropy in the same 

wavenumber range. We found that a plateau in the compensated three-dimensional energy spectrum at 

kr) & 0.1 0.2, commonly used to infer the Kolmogorov constant from the compensated three-dimensional 

energy spectrum, actually does not represent proper inertial range behavior. Rather, a proper, if still 
approximate, inertial range emerges at kr] « 0.02 — 0.05 when R\ increases beyond 140. The new simulations 
indicate proportionality constants C\ and C in the one- and three-dimensional energy spectra respectively 
about 0.60 and 1.62. If the turbulence were perfectly isotropic then use of isotropy relations in wavenumber 
space (Ci = 18/55 C) would imply that Ci « 0.53 for C = 1.62, in excellent agreement with experiments. 
However the one- and three-dimensional estimates are not fully consistent, because of departures (due to 
numerical and statistical limitations) from isotropy of the computed spectra at low wavenumbers. The 
inertial scaling of structure functions in physical space is briefly addressed. Since DNS is still restricted to 
moderate Reynolds numbers, an accurate evaluation of the Kolmogorov constant is very difficult. We focus 
on providing new insights on the interpretation of Kolmogorov 1941 similarity in the DNS literature and do 
not consider issues pertaining to the refined similarity hypotheses of Kolmogorov (K62). 

Key words. Kolmogorov constant, energy spectrum, numerical simulations 

Subject classification. Fluid Mechanics 

1. Introduction. Inertial-range behavior as postulated by Kolmogorov’s (K41) similarity hypotheses 
[1] is widely regarded as a fundamental characteristic of turbulence at high Reynolds number. In particular, in 
the inertial range of intermediate scales the K41 result for the one- dimensional longitudinal energy spectrum 
is given by 

(1) E u (ki) = C 1 {ef /3 k; 5/3 , 

where k\ is the longitudinal wavenumber, C\ is known as the Kolmogorov constant, and (e) is the mean 
dissipation rate. Although K41 theory is subject to intermittency corrections associated with dissipation rate 
fluctuations, such effects are primarily manifested in higher-order statistics. Indeed, intermittency effects on 
the second-order energy spectrum exponent are believed to be small and hardly measurable (Kolmogorov 

[2] , Frisch [3]), while at the same time may contribute to a persistence Reynolds number dependence for 
higher-order structure functions (L’vov & Procaccia [4]). 
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The classical view of the “five-thirds” scaling law above, with substantial experimental support (e.g., 
Monin &; Yaglom [5], Sec. 23.3), is that C\ has an universal value at asymptotically high Reynolds number. 
Recently, however, there is renewed debate on the universality of Ci, in part because of new measurements 
at high Reynolds number (Praskovsky & Oncley [6]) and a subsequent new similarity theory (Barenblatt & 
Goldenfeld [7]) that suggested a persistent Reynolds number dependence even at high Reynolds numbers. 
On the other hand, the conclusion from a new and very extensive survey of experimental data by Sreeni- 
vasan [8] is that, taken collectively, measurements do not support such a dependence for the Kolmogorov 
constant. Sreenivasan found that the value of C\ averaged over the many different experiments cited is about 
0.53 ± 0.055, although some corrections for the estimation of dissipation using local isotropy assumptions in 
experiments may be warranted. It is still possible (Barenblatt & Goldenfeld [7], Sreenivasan [8]), though, 
that controlled experiments for a single geometry over a wide range of Reynolds numbers may reveal signifi- 
cant trends otherwise masked by experimental scatter; new measurements of this nature (suggesting strong 
Reynolds number dependence) have been reported by Mydlarski & Warhaft [9], 

This paper is motivated by the survey of experimental data noted above, and will focus on similar 
issues arising in numerical simulations of isotropic turbulence. We first review the basis for estimating the 
Kolmogorov constant from previous studies, and then present new results from direct numerical simulations 
(DNS) over a range of Reynolds numbers. Because of Reynolds number considerations, an accurate evaluation 
of the Kolmogorov constant by DNS is admittedly very difficult. As such, we limit ourselves to the task of 
providing new insights on the interpretation of K41 similarity in the DNS literature. For similar reasons, 
we shall not consider the use of DNS to study issues pertaining to the refined similarity hypotheses of 
Kolmogorov (K62) [2] (see, e.g., Chen et al [10], Wang et al [11] and others). 

In numerical simulations, K41 similarity is frequently discussed in terms of the three-dimensional energy 
spectrum function E(k) (where k is the wavenumber magnitude). If the turbulence at scale size 1 jk is 
isotropic then a kinematic constraint relating one- and three-dimensional spectra is 


( 2 ) 



d_ 

dk 


(ldE xl (k) 
\k dk 


) 


Substitution of Eq. 1 into Eq. 2 implies that in the inertial range 


(3) 


E(k) = C(ef /3 k~^ 3 , 


where C = 55/18 C\. In principle, therefore, one may obtain either C or C\ (from three- and one-dimensional 
spectra respectively) and infer the other using isotropy relations. Values of the Kolmogorov constant C 
(estimated using either of these approaches) cited in a number of recent high resolution numerical studies 
by other authors [11-21] are collected in Table I. Some relevant results from large-eddy simulations (LES) 
are included also. (Note that we have quoted values of C from the highest-resolution data in each of these 
references, since from the viewpoint of present paper we should not infer the Kolmogorov constant from 
low-resolution simulations, especially in the older data.) 

It may be noticed that most of the values displayed in Table I are higher than the value 1.619, which is 
55/18 of the experimental average of 0.53 for C\. Note that, however, Zhou [14] found that the Kolmogorov 
constant C = 1.5 from calculations of the spectra energy flux based on the self- similarity condition of the 
triadic energy transfer function, in the form T(k,p, q) = a z T(ak,ap,aq) (where a is a scaling constant). A 
commonly used procedure (e.g. Kerr [17], Jimenez et al [20]) for estimating the value of C is to plot the 
“compensated” three-dimensional energy spectrum 

(4) ip(k) = E(k)(ey 2/3 k 5 ' 3 
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as a function of wavenumber normalized by the Kolmogorov scale ( 77 ), and to interpret C as the height of a 
“plateau”. However, it should be noted that whereas a flat region in the compensated spectrum implies k~ 5 ^ 3 
behavior, the observation of a k~ 5//3 scaling range is not by itself a sufficient condition for an inertial range. 
It is important that isotropy also be attained in the wavenumber range which displays k ~ 5//3 behavior. (If 
the isotropy requirement were relaxed then in principle one would have a different value of the “Kolmogorov 
constant” in each direction, and the concept of universality would be lost.) Indeed, as will be seen in the 
rest of this paper, our new results demonstrate that deviations from isotropy can contribute to values of C 
that appear to be too high. 

The major limitation in using DNS to examine inertial-range dynamics is, of course, the difficulty in 
attaining high Reynolds numbers. However, recent advances in massively parallel computing have shown 
significant promise. Our simulations were performed using a parallel implementation [ 22 ] of the well-known 
Fourier pseudo-spectral algorithm of Rogallo [23] on the IBM SP at the Cornell Theory Center. The highest 
grid resolution used is 512 3 , with a Taylor-scale Reynolds number {R\) about 240 averaged over about 
four eddy- turnover times. Whereas this Reynolds number is not high compared to some recent laboratory 
experiments (R\ 473 in Mydlarski & War haft [9]), it is about the same as (or slightly higher than) the 
highest values reported in the DNS literature (for example, R\ 218 in Cao et ai [24]). Numerical results 
on spectra as well as structure functions (to which Kolmogorov [1] originally referred) are given in the next 
section. In order to characterize issues of Reynolds number dependence and attempt to quantify a minimum 
Reynolds number threshold above which inertial-range behavior could be expected, we present data at five 
different grid resolutions. 

To maintain a stationary state so that results may be averaged over time at the highest Reynolds number 
possible using a given number of grid points, it is usual to apply numerical forcing at the large scales. In the 
literature this has been done in several different ways, such as adding a stochastic forcing term (Eswaran 
Sz Pope [25]), holding the energy in the lowest- wavenumber shells fixed while allowing phase information 
to evolve (Chen et ai [10], Sullivan et ai [26]), or by introducing a negative viscosity in these shells 
(Jimenez et ai [20]). Alternatively, LES [12-15], with and without forcing, as well as simulations with a 
hyperviscosity [16] can also be used to achieve higher Reynolds numbers. However, it is generally believed 
that (e.g., see Ref. 20 ), the precise manner of forcing- provided it is applied to the largest scales in the 
flow— has no systematic effects on the inertial- range energy spectrum, nor on the statistical character of 
the small scales. Whereas in this work we have used the scheme of Eswaran & Pope [25], we have checked 
that similar calculations in which the energy in the first couple of wavenumber shells is fixed do not warrant 
different conclusions. 

2. Results. We present spectra and structure functions from simulations of forced stationary isotropic 
turbulence averaged over relatively long time periods. Major simulation parameters are summarized in 
Table II. For the less expensive 64 3 and 128 3 runs the averaging time (T) shown is the aggregate of multiple 
simulations of shorter duration. In all cases the small scales are considered to be well resolved, as measured 
by k max r] « 1.5, where k max is the highest wavenumber represented by the grid points. The non-dimensional 
quantity (e)Li/u f3 (where L\ is the longitudinal integral length scale derived from the longitudinal velocity 
correlation, and u f is the r.m.s. velocity) in the last column of this table is of interest in the scaling of (c) 
with u f3 /L\ using energy cascade arguments; it appears to approach a constant at high Reynolds numbers. 
This trend is similar to that found in the simulations of Jimenez et al. [20], Wang et ai [ 11 ] and Cao et 
ai [24], and to that in experiments (Sreenivasan [27], Fig. 1 ). The differences among different “asymptotic” 
numerical values are due in part to differences in the definition of integral scales and in flow conditions in 
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experiments and numerical simulations. 

Figure 1 shows the three-dimensional compensated spectrum ip(k) in Kolmogorov scaling. (Note: in 
these and all other figures fines A — E refer to data at five different grid resolutions from 64 3 to 512 3 
respectively, as fisted in Table II.) Except for small spectral turn-ups at the high wavenumber end resulting 
from residual aliasing errors, small-scale universality is unambiguously achieved, even at the lowest Reynolds 
number ( R\ 38 in 64 3 simulations). Two relatively flat (Le. A; -5 / 3 ) regimes of limited extent can be seen: 
namely at A:77 ^0.1 — 0.2 in the form of a bumpy “plateau” which is seen at all grid resolutions, and in the 
low-wavenumber range krj ~ 0.02 — 0.05 which is captured only at higher (256 3 and beyond) resolutions. 
If the former were taken to represent the inertial range then one would obtain a Kolmogorov constant C 
greater than 2.0, as reported by a number of other authors (see Table I). On the other hand, it is clear that 
the level of 'ip(k) in the range kr} « 0.02 — 0.05 agrees well with experimental data, if one uses the isotropy 
relation C = 55/18 C\ , which would imply C\ = 18/55 C & 0.53. It is our objective in the analyses below 
to establish that this lower-wavenumber region, rather than the plateau, represents (the beginnings of) a 
proper inertial range. 

In order that inertial range dynamics can be independent of both the large-scale energetics and viscous 
dissipation, a wide scale separation must exist between the peaks of the energy and dissipation spectra in 
wavenumber space. The peak of E(k) occurs in the lowest two wavenumber shells (krj % 0.01 for the 512 3 
data), whereas from a plot of the dissipation spectrum we find that the peak of D(k) = 2 vk 2 E(h) occurs 
at krj « 0.17. This location of the dissipation spectral peak is virtually the same as found by Wang et 
ai [11]. Since this almost coincides with the “plateau” in Fig. 1 it is clear that the plateau should not be 
taken as an indication of inertial- range behavior. In fact, this plateau may be identified with the so-called 
“bottleneck” phenomenon that has been discussed in experimental (Saddoughi & Veeravalli [28]), theoretical 
(Falkovich [29]), and numerical (Borue & Orszag [30]) work. According to this theory, viscous effects on 
spectral transfer causes an increase (relative to k~ 5//3 in the energy spectrum at intermediate wavenumbers. 
Nonlocal interactions among Fourier modes well separated in scale (Brasseur & Wei [31]) are also believed 
to play an important role (Herring et ai [32]). We also note that Brasseur & Wei argued that there should 
be about a decade in scale separation between the wavenumber where the inertial range “ends” and the 
dissipation peak, which is more than that (about 3 to 4) found in DNS (Wang et ai [11] and this paper). 
A full understanding of this latter issue awaits future high-resolution simulations and experimental data. 

A more direct comparison with experiment may be made by considering the compensated longitudinal 
energy spectrum, i.e. ipi(ki) = Ei\(ki)(e)~ 2 ^k\^ , which is shown in Fig. 2. Because this spectrum drops 
off very rapidly at high wavenumbers, we have used log- linear scales to highlight the function values at any 
wavenumber ranges where Vfi is approximately constant. Two such ranges can be seen, corresponding to 
those for if) but in each case occurring at lower wavenumbers. The observation that the peak of 
occurs at k\rj = 0.05 is similar to high Reynolds number measurements in both boundary layers (Saddoughi 
& Veeravalli [28]) and grid turbulence (Mydlarski & Warhaft [9]). It is also apparent that the value of 
C\ inferred from this figure is about 0.60, which is somewhat higher than, but still relatively close to the 
experimental average [5] of 0.53. 

From Figs. 1 and 2 we may conclude that the best values obtained for C\ from the three- and one- 
dimensional spectra are about 0.53 and 0.60 respectively. These values (especially the former) are closer to 
experiment than those cited previously in the literature. Yet there is evidently some inconsistency between 
these two estimates for C\. These differences are due to deviations from isotropy in wavenumber space, as 
further studied below. 
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In addition to En(k\) we have also computed the transverse energy spectrum, for which the classical 
inertial- range result is 

(5) E 2 2 (h) = C[(e) 2/3 k; 5/3 , 

where isotropy requires C\ = |(7i . The compensated transverse spectrum (Fig. 3) exhibits a similar but 
stronger bump at k\T] ~ 0.07. We also find that the ratio £^ 2(^1 )/£n(^i) Is close to 4/3 in the range 
kit) & 0.02 — 0.04, but increases steadily with wavenumber beyond this range. 

If a five-thirds scaling holds for both one- and three-dimensional spectra over the same wavenumber 
range then the functions E(k) and En(k) should be proportional within this range. Figure 4 shows the 
spectral ratio E(k)/Eu(k), compared with the classical inertial range value of 55/18. It may be seen that, 
at the wavenumber range (kr) ~ 0.1 — 0.2) corresponding to the plateau in Fig. 1, this ratio is considerably 
higher than 55/18 and in fact increases roughly in proportion to kr]. This observation provides further 
evidence that the plateau in ^ does not represent an inertial range. The transition from behavior without 
a level region to one at a value close to 55/18 (at R\ 240 at 512 3 resolution) occurs between R\ 90 (on an 
128 3 grid) and R\ 140 (for 256 3 ). The ratio C/Ci inferred in this manner is seen to be somewhat less than 
55/18. 

It should be noted that isotropy of the spectral tensor is a key requirement for the relation C = 55/18Ci. 
Although the simulations are of (nominally) isotropic turbulence, departures of the computed spectra from 
isotropy are, in fact, not totally unexpected. It is well-known that sampling limitations arise in the lowest 
few wavenumber shells because relatively few samples of the large scales are present in a solution domain 
of finite size. Furthermore, since the solution domain in physical space is a cube (rather than a sphere), 
the averaged statistics are at best invariant among the three Cartesian coordinate axes, but do not satisfy 
the stricter requirement of no preferential orientation in three-dimensional space. As an example we may 
note that whereas the periodicity length of the flow along the coordinate axes is equal to the length of each 
side of the solution domain, it is longer if measured along directions inclined to the axes. This effect is felt 
primarily in spatial correlations over large separations in space, which correspond to low wavenumbers in 
Fourier space. 

In Figs. 6 and 7 these structure functions are shown in Kolmogorov scaling in order to compare with 
the theoretical proportionality constants C 2 ~ 4.02 C\ and —4/5. It is clear that the highest Reynolds 
number data shown agree well with classical inertial range results. Numerical values of — {t) r 

intermediate resolutions also fit in well within the Reynolds number trend suggested by new measurements 
in grid turbulence over a range of Reynolds numbers (Sreenivasan & Dhruva [35]). 

It is perhaps worth noting that the highest-rcsolution (51 2 3 ) simulation reported is at a higher Reynolds 
number and are averaged over a greater number of large-eddy turnover times than several other studies 
reporting 512 3 results (e.g., Jimenez et al [20], Chen et al [10], Wang et al [11]) for stationary 
isotropic turbulence. The spectra and structure functions obtained from these and the present simulations 
demonstrate that, with the latest advances in massively parallel computing, issues concerning inertial- range 
similarity in DNS can now be addressed in a more reliable manner than possible before. 

3. Conclusions . We have presented new results on the Kolmogorov scaling of energy spectra and 
structure functions in the inertial range, from direct numerical simulations of stationary isotropic turbulence 
ranging from R\ 38 (on a 64 3 grid) to R\ 240 (on 512 3 ). It is pointed out that a plateau in the compen- 
sated three-dimensional energy spectrum at kr} ~ 0. 1-0.2 commonly used to infer the Kolmogorov constant 
from the compensated three-dimensional energy spectrum in fact does not represent proper inertial range 
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behavior. Instead, a proper (if still approximate) inertial range emerges at hr) ~ 0.02 — 0.05 when R\ in- 
creases beyond 140. We find that the proportionality constants C and C\ in the three- and one-dimensional 
compensated energy spectra are about 1.62 and 0.60 respectively. These values are closer to experimental 
data than reported in most previous numerical simulations. In particular, if isotropy relations are used then 
we may infer from the three-dimensional spectra that C\ — 18/55 C ~ 0.53, in excellent agreement with 
experimental data (Sreenivasan [8]). However, the ratio C/C\ in our results differs from the theoretical value 
of 55/18, because of significant departures from isotropy in the computed spectra at the wavenumber range 
where inertial- range behavior is otherwise reasonably well approximated. Results on second- and third-order 
structure functions over a range of Reynolds numbers further suggest that the simulation database that we 
have accumulated should be useful for investigating other aspects of similarity scaling and Reynolds number 
dependence. 

We emphasize that a scaling is in itself not a sufficient indicator of inertial- range behavior. To 

achieve a strictly isotropic inertial range in direct numerical simulations requires that the inertial scales be 
small compared to the size of the solution domain. Whereas it is difficult to meet this requirement well, it 
seems clear that high-resolution simulations using the techniques of massively parallel computing are very 
helpful. 
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Fig. 1. Compensated three-dimensional energy spectrum tp(k) versus Kolmogorov- scaled wavenumber hr), at five different 
Reynolds numbers and grid resolutions (see Table II). Dashed horizontal lines at levels 1.619 = (55/18)0.53 and 2.0 are drawn 
for reference. 
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Table 1 


Kolmogorov i 

constant computed from DNS and LES a 


authors 

method 

Grid 

R\ 

C (ID) 6 

C (3D) 

LesieurSzRogallo [12] 

Forced LES C 

128 3 



1.5-1. 8 

Chasnov [13] 

Forced LES d 

128 3 



2.1 

Zhou [14] 

LES e 

256 3 

- 

- 

1.5 

SheSzJackson [15] 

LES' 

128 3 

-- 


1.88 

BaruehOrszag [16] 

LES 9 

256 3 

- 

1.2-1. 7 

1.2 -1.4 

Kerr[ 17] 

Forced DNS 

128 3 

82 


2 

Vincent&cM eneguzzi[ 18] 

Forced DNS 

240 3 

150 


2 

Sanda[ 19] 

Forced DNS 

256 3 

120 


2 

Jimenez et al [20] 

Forced DNS 

512 3 

170 - 

0.91 1.22 

2. 

Wang et aL[ll] 

Forced DNS* 

512 3 

190 

1.68 

1.5-2. 

Hosokawa et al. [21] 

DNS 

512 3 

160 

2.1 

- 

Present study 

Forced DNS 

512 3 

240 

1.83 

1.62 


Footnotes 

a We included most of the recent numerical data on the Kolmogorov constant. 

b C (ID) is obtained as (55/18)Ci, assuming isotropy in the inertial range. (Hosokawa et al. used 0.76C 2 , 
from the structure function.) 
c Spectral LES with a spectral eddy viscosity. 
d Spectral LES with an eddy viscosity and a stochastic force. 

e Zhou used simulation databases from a constrained energy simulation where a “five-thirds” spectrum 
is maintained. Furthermore, the self-similarity condition is used to compute an ideal Kolmogorov energy 
transfer function. The Kolmogorov constant is determined from the flux. 

' She and Jackson used the constrained Euler system simulation. The Kolmogorov constant is determined 
from the flux. 

9 Borue and Orszag reported their decaying DNS at 256 3 resolution along with a hyperviscosity. Both C 
(ID) and C (3D) are estimated by the authors of present paper from Fig. 4 of Ref. 16. The contributions 
from the “bumps” are ignored. 

h Wang et al argued that the Kolmogorov constant is more accurate when estimated from C (ID) followed 
by using the isotropy relation. 
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Table 2 

Major simulation parameters in the new results. 


Grid 

Rx 

T/T e 

{t)L\/u l 

64 3 

38 

146 

0.693 

128 3 

90 

139 

0.499 

256 3 

140 

9.8 

0.475 

384 3 

180 

3.2 

0.419 

512 3 

240 

3.9 

0.416 


Note: T is the total averaging time period, and Te is the eddy-turnover time (the ratio of the longitudinal 
integral length scale, L\ y to the r.m.s. velocity fluctuation, u') 


10 




Fig. 2. Compensated one- dimensional energy spectrum versus Kolmogorov- scaled wavenumber kirj, at five different 
Reynolds numbers and grid resolutions (see Table II). The reference lines are drawn at levels 0.53 and 0.60. 



Fig. 3. Same as Fig. 2, but for the transverse spectrum. The reference lines are drawn at levels 0.707 = (4/3)0.53 and 0.8. 
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Fig. 4. The ratio E(k)/E\i(k) versus Kolmogorov- scaled wavenumber. The horizontal line is at 55/18, whereas the sloping 
line has slope unity. Because this ratio attains very large values at high wavenumbers , only the lower half of the wavenumber 
range in each simulation is shown. 



Fig. 5. The isotropy coefficient I(k\) (see Eq. 6) versus Kolmogorov- scaled one- dimensional wavenumber. (A value of 
unity signifies spectral isotropy.) 
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Fig. 6. Kolmogorov scaling of the second- order longitudinal structure function versus r/eta. The horizontal line is at 
2.13 = (4.02)(0.53); the sloping dashed line indicates universal behavior of the small scales, as (l/15)(r/r7) 4/3 . 



Fig. 7. Kolmogorov scaling of the third order longitudinal structure function versus r/eta. The horizontal line is at 0.8. 
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